%\documentclass{ctexart} 使用 ctex 字体
\documentclass{ctexart} 

\usepackage{xltxtra} % 让 xelatex 支持UTF
\usepackage{graphicx} % 图形库
\usepackage{amsmath} % AMS 的数学符号
\usepackage{pgf,pgfarrows,pgfnodes}
\usepackage{listings}

\title{step3的描述}


\author{ 梁越 \\ 信息与计算科学 3220100159}

\begin{document}


\maketitle

\section{step3}
\subsection{有限元方法的基本设置}
这是我们实际使用有限元计算的第一个例子。我们将求解边界值为零但右侧非零的简单泊松方程：
\begin{equation}
-\mathrm{d}u=f in\Omega ,
\end{equation}
\begin{equation}
u=0 on\partial \Omega 
\end{equation}
我们将在单位平方上解决这个方程，$\Omega=[0,1]^{2} $，为此你已经学习了如何生成一个网格step-1和step-2。在这个程序中，我们也将只考虑特殊情况f(x)=1，并在下一个教程程序步骤4中回到如何实现更一般的情况。
如果你学过有限元法的基础知识，你会记得我们需要采取的步骤，用有限维近似来近似解。具体来说，我们首先需要推导出上面方程的弱形式，我们可以通过从左边乘以一个测试函数$\phi $来得到它(我们将回到从左边乘而不是从右边乘的原因)，并对$\Omega$进行积分:
\begin{equation}
-\int_{\Omega}^{} \varphi \mathrm{d}u = \int_{\Omega}^{} \varphi f.
\end{equation}
这可以通过以下部分进行集成：
\begin{equation}
\int_{\Omega} \nabla\phi \cdot \nabla u - \int_{\partial\Omega} \phi n \cdot \nabla u = \int_{\Omega} \phi f.
\end{equation}
在这里，我们使用了常见的符号表示：$(a, b) = \int_{\Omega} ab $。该问题要求找到一个函数 (u) ，对于所有来自适当空间（在这里为$(H^1)$空间）的测试函数$\varphi$ ，使得下述陈述成立：
\begin{equation}
\int_{\Omega} \nabla \varphi \cdot \nabla u - \int_{\partial \Omega} \varphi n \cdot \nabla u = \int_{\Omega} \varphi f
\end{equation}
在一般情况下，我们无法在计算机上找到这样的函数。因此，我们寻求一个近似解 $u_h(x) = \sum_j U_j \varphi_j(x)$，其中 $(U_j)$ 是待确定的未知系数（问题的“自由度”），$\varphi_i(x)$ 是我们将使用的有限元形函数。为了定义这些形函数，我们需要以下内容：
    用于定义形函数的网格。在第一步和第二步中，您已经了解了如何生成和操作描述网格的对象。
    描述我们希望在参考单元上使用的形函数的有限元。在第二步中，我们已经使用了类型为 $FE_Q<2>$ 的对象，它代表了通过插值支持点定义形函数的通常的拉格朗日元素。最简单的情况是 $FE_Q<2>(1)$，它使用一次多项式。在二维中，它们通常被称为双线性函数，因为它们在参考单元的两个坐标中都是线性的。（在一维中，它们将是线性的，在三维中将是三线性的；然而，在deal.II的文档中，我们通常不做这个区分，而只是称这些函数为“线性”。）
    ADoFHandler 对象，该对象枚举了网格上的所有自由度，以有限元对象提供的参考单元描述为基础。在第二步中，您已经了解了如何执行此操作。
    一个映射，用于告诉我们如何从参考单元上有限元类定义的形函数获取实际单元上的形函数。默认情况下，除非您明确指定，deal.II将使用（双线性、三线性）映射进行此操作，因此在大多数情况下，您不必担心这一步。
通过上述步骤，我们现在有了一组函数 $\varphi_i$，我们可以定义离散问题的弱形式：找到一个函数 $u_h$ ，即找到上述提到的展开系数 $U_j$，使得
\begin{equation}
(\nabla \varphi_i, \nabla u_h) = (\varphi_i, f), \quad i=0,\dots,N-1
\end{equation}
请注意，这里我们遵循从零开始计数的惯例，与C和C++中常见的方式相符。如果插入表示 $u_h(x) = \sum_j U_j \varphi_j(x)$ ，则可以将该方程重写为线性系统，然后观察到
\begin{equation}
(\nabla \varphi_i, \nabla u_h) = (\nabla \varphi_i, \nabla \left[\sum_j U_j \varphi_j\right]) = \sum_j (\nabla \varphi_i, \nabla [U_j \varphi_j]) = \sum_j (\nabla \varphi_i, \nabla \varphi_j) U_j
\end{equation}
这样，问题可以表示为：找到一个向量 $\mathbf{U}$ 使得
\begin{equation}
\mathbf{AU} = \mathbf{F},
\end{equation}
其中矩阵 $\mathbf{A}$ 和右手边向量 $\mathbf{F}$ 定义为
\begin{equation}
A_{ij} = (\nabla \varphi_i, \nabla \varphi_j), \quad F_i = (\varphi_i, f).
\end{equation}
\subsection{我们应该在左边还是右边乘以一个测试函数？}
在我们继续描述这些量如何计算之前，需要注意如果我们将原始方程从右边而不是左边乘以一个测试函数，则会得到一个形式为
\begin{equation}
(\mathbf{U}^T \mathbf{A} = \mathbf{F}^T)
\end{equation}
的线性系统，其中 $\mathbf{F}^T$ 是行向量。通过转置这个系统，这当然等价于求解
\begin{equation}
(\mathbf{A}^T \mathbf{U} = \mathbf{F})
\end{equation}
在这里，由于 $\mathbf{A} = \mathbf{A}^T$，所以与上面的情况相同。但是一般来说是不同的，在避免任何混淆的情况下，经验表明习惯上将方程从左边而不是右边乘以一个测试函数（这在数学文献中经常这样做）可以避免一类常见的错误，因为矩阵会自动正确，不需要在比较理论和实现时转置。请参阅这个教程中的第一个示例（第9步），其中我们有一个非对称的双线性形式，从右边或左边乘以测试函数会有所不同。
\subsection{计算矩阵和右手边向量}
现在我们知道我们需要什么（即：存储矩阵和向量的对象，以及计算 $A_{ij}$ 和 $F_i$ 的方法），我们可以看一下如何实现：

$\mathbf{A}$、$\mathbf{U}$ 和 $\mathbf{F}$ 的对象是 SparseMatrix 和 Vector 类型的，我们将在下面的程序中看到用于解线性系统的哪些类。
我们需要一种方法来进行积分。在有限元方法中，最常用的方法是使用数值积分，即用每个单元上一组点的加权和来代替积分。也就是说，我们首先将对 $\Omega$ 的积分分解为所有单元上的积分，
\begin{equation}
A_{ij} = (\nabla \varphi_i, \nabla \varphi_j) = \sum_{K \in \mathcal{T}} \int_K \nabla \varphi_i \cdot \nabla \varphi_j, \quad F_i = (\varphi_i, f) = \sum_{K \in \mathcal{T}} \int_K \varphi_i f,
\end{equation}
然后通过数值积分近似每个单元的贡献：
\begin{equation}
A_{Kij}F_{Ki} = \int_K \nabla \varphi_i \cdot \nabla \varphi_j \approx \sum_q \nabla \varphi_i(x_{Kq}) \cdot \nabla \varphi_j(x_{Kq}) w_{Kq}, \quad F_{Ki} \approx \sum_q \varphi_i(x_{Kq}) f(x_{Kq}) w_{Kq},
\end{equation}
找到向量 $\mathbf{U}$，使得
\begin{equation}
\mathbf{AU} = \mathbf{F},
\end{equation}
其中矩阵 $\mathbf{A}$ 和右侧向量 $\mathbf{F}$ 定义为
\begin{equation}
A_{ij} = (\nabla \varphi_i, \nabla \varphi_j), \quad F_i = (\varphi_i, f).
\end{equation}
为了计算这些量，我们需要以下步骤：

    确定积分点的位置 $x_{Kq}$ 和权重 $w_{Kq}$，通常通过参考单元上的积分点和权重与映射函数（Mapping类）相结合得到。可以使用Quadrature基类来描述参考单元上积分点和权重，而QGauss类提供了高斯积分的实现。

    使用FEValues类来计算在单元K上的形函数值 $\varphi_i(x_{Kq})$。FEValues类需要有限元对象来描述参考单元上的形函数，以及映射对象和积分点的信息。FEValues类提供了在积分点上计算形函数值和导数等所需的信息。

    FEValues类是装配过程中的关键类，它将有限元、积分点和映射对象相结合，提供在特定单元K上所需的有限信息集。在装配线性系统时，你将看到FEValues类的实际应用。

最后，通过迭代求解器解决线性系统后，使用DataOut类创建输出文件，并可以使用常见的可视化程序进行可视化。
\footnote{前面对任何有限元实现的所有重要步骤的概述在交易中都有对应的部分。II：该库可以自然地分为多个“模块”，涵盖刚刚概述的基本概念。您可以通过此页面顶部的选项卡访问这些模块。最基本的概念组概述也可在交易.II手册首页 .}
\subsection{关于实施}
虽然这是使用有限元方法可以求解的最简单的可能方程，但该程序显示了大多数有限元程序的基本结构，并且还充当了几乎所有以下程序基本遵循的模板。具体来说，该程序的主类如下所示：
\begin{lstlisting}[language=C++]
classStep3
{
public:
Step3 ();
voidrun();
private:
voidmake_grid ();
voidsetup_system ();
voidassemble_system ();
voidsolve ();
voidoutput_results ()const;
Triangulation<2>triangulation;
FE_Q<2>fe;
DoFHandler<2>dof_handler;
SparsityPatternsparsity_pattern;
SparseMatrix<double>system_matrix;
Vector<double>solution;
Vector<double>system_rhs;
};
\end{lstlisting}
这遵循了面向对象编程的原则数据封装也就是说，我们尽力将这个类的几乎所有内部细节隐藏在外部无法访问的私有成员中。

让我们从成员变量开始：这些变量遵循我们在上述要点中概述的构建块，即我们需要一个三角测量和一个DORHANDLER对象，以及描述我们要使用的形状函数类型的有限元对象。第二组对象与线性代数有关：系统矩阵、右手边以及解向量，以及描述矩阵稀疏模式的对象。这就是该类所需的全部内容（以及任何静态PDE求解器所需的基本要素），并且需要在整个程序中生存下来。与此相反 FE值装配所需的对象仅在整个装配过程中是必需的，因此我们在执行此操作的函数中将其创建为局部对象，并在其末尾再次销毁它。

其次，让我们看看成员函数。这些也已经形成了几乎所有以下教程程序都将使用的通用结构：
\verb|make_grid()|：这就是我们可以称之为预处理函数顾名思义，它设置存储三角剖分的对象。在后面的示例中，它还可以处理边界条件、几何图形等。
\verb|setup_system()|：这就是解决问题所需的所有其他数据结构的设置函数。特别是，它将初始化DOF HANDLER，并正确调整与线性代数有关的各种对象的大小。此函数通常与上面的预处理函数分离，因为在依赖时间的程序中，每当网格自适应细化时，至少每隔几个时间点调用一次此函数（我们将在步骤-6). 另一方面，在上述预处理函数中设置网格本身只在程序开始时进行一次，因此，它被分为自己的函数。
\verb|assemble_system()|：这就是计算矩阵内容和右侧内容的地方，如上文引言中详细讨论的那样。由于对这个线性系统进行操作在概念上与计算其条目非常不同，因此我们将其与以下函数分开。
\verb|solve()|:这就是我们计算线性系统$ AU = F $的解的函数。在当前的程序中，这是一个简单的任务，因为矩阵是如此简单，但是当问题不再那么简单时，它将成为程序大小的重要组成部分(例如，一旦您对库有了更多的了解，请参阅步骤20，步骤22或步骤31)。
\verb|output_results()|：最后，当你计算出一个解时，你可能想用它做点什么。例如，你可能想要以一种可以可视化的格式输出它，或者你可能想要计算你感兴趣的量：例如，热交换器中的热通量，机翼的空气摩擦系数，最大桥梁荷载，或者只是数值解在某一点的值。因此，此函数是对解决方案进行后处理的位置。
所有这些都由单个公共函数（而不是构造函数）结合在一起，即\verb|run()|功能。它是从创建此类型对象的位置调用的，也是按正确顺序调用所有其他函数的。将此操作封装到运行函数，而不是从中调用所有其他函数\verb|main()|确保可以更改此类中关注点分离的实现方式。例如，如果其中一个函数变得太大，您可以将其拆分为两个，因此您唯一需要关心更改的地方是在这个类中，而不是在其他任何地方。

如上所述，您将在下面的许多教程程序中再次看到这种通用结构——有时在函数名称的拼写中会有变体，但本质上是按功能分离的顺序。
\subsection{关于类型的注释}
deal.II通过命名空间类型中的别名定义了许多整数类型。特别是，在这个程序中，您将看到类型：\verb|global_dof_index|在几个地方：整数类型，用于表示自由度的全局索引，即在三角测量上定义的dofhandlerobject内特定自由度的索引(与特定单元格内特定自由度的索引相反)。对于当前的程序(以及几乎所有的教程程序)，您将在全球范围内拥有几千到几百万个未知量(并且，对于q1元素，您将在2d中每个单元格上有4个，在3d中有8个)。因此，允许为全局DoF索引存储足够大的数字的数据类型是\verb|unsigned int|考虑到它允许存储0到略多于40亿的数字（在大多数系统上，整数是32位的）。事实上，这就是类型\verb|global_dof_index|

那么，为什么不使用\verb|unsigned int|？在7.3版之前，我一直这样做。然而，deal.II支持非常大的计算（通过中讨论的框架步骤40)当分布在数千个处理器上时，可能有超过40亿个未知项。因此，在某些情况下\verb|unsigned int|不够大，我们需要一个64位无符号整数类型。为了实现这一点，我们引入了类型\verb|global_dof_index|默认情况下，其定义如下\verb|unsigned int|然而，可以将其定义为无符号长长整型如有必要，在配置期间传递特定标志（请参阅自述文件）。

这包括技术方面。但也有一个文档目的：如果你看到一个地方使用数据类型，那么库中的任何地方以及构建在其上的代码类型：\verb|global_dof_index|，您立即知道所引用的数量实际上是一个全局自由度指数。如果我们只是使用了\verb|unsigned int|（也可以是局部索引、边界指示器、材质id等）。立即了解变量所指的内容也有助于避免错误：很明显，如果看到类型为的对象，一定会有错误类型：：\verb|global_dof_index|分配给类型为的变量类型：：\verb|subdomain_id|，即使它们都由无符号整数表示，编译器也不会因此而抱怨。
在更实际的术语中，这种类型的存在意味着在组装期间，我们创建a4×4matrix(在2d中，使用aq1元素)我们当前所处的单元格的贡献，然后我们需要将该矩阵的元素添加到全局(系统)矩阵的适当元素中。为此，我们需要获得当前单元局部自由度的全局指标，为此我们将始终使用以下代码片段:
\begin{lstlisting}[language=C++]
cell->get_dof_indices (local_dof_indices);
\end{lstlisting}
此处\verb|local_dof_indices|声明为
\begin{lstlisting}[language=C++]
std::vector<types::global_dof_index> local_dof_indices $fe.dofs_per_cell$;
\end{lstlisting}
这个变量的名称可能有点用词不当，它代表“当前单元格中本地定义的那些自由度的全局索引”，但保存此信息的变量在整个库中都是这样命名的。\footnote{类型：$global_dof_index$不是此命名空间中定义的唯一类型。相反，有一个完整的家庭，包括类型：$subdomain_id$ ,类型：$boundary_id$和类型：$material_id$所有这些都是整型数据类型的别名，但如上所述，它们在整个库中使用，以便更容易识别变量的意图，如有必要，可将实际类型更改为更大的类型}
\section{step-3问题 参数 图片}

\subsection{简介}
本报告中的工作基于deal.II库的示例程序step-3。该程序主要用于解决二维Poisson方程的数值计算问题。具体而言，我们将通过离散化有界区域Ω上的Poisson方程来计算解。在边界上，我们会施加适当的边界条件，以得到解的准确解析表达式
\subsection{数学模型}
首先，让我们定义$\Omega $为一个单位正方形$(0,1)x(0,1)$。Poisson方程的数学模型如下所示
\begin{equation}
        \Delta u = f \text{ in } \Omega
\end{equation}
\begin{equation}
        u = g \text{ on } \partial \Omega
\end{equation}
其中，$\Delta$是Laplace算子，u是待解函数，f是源项函数，g是边界条件。在本示例中，我们将考虑以下参数设置：
\begin{equation}
f(x,y) = 1
\end{equation}
\begin{equation}
g(x,y) = x*(1-x)*y*(1-y)
\end{equation}
\subsection{数值分析方法}

我们使用有限元方法对Poisson方程进行数值求解。首先，我们将Ω分割成一组小的三角形（网格）。然后，在每个三角形上构建有限元空间，其中解u在每个元素中用一组基函数展开。最后，通过求解线性代数方程组来计算解的系数。


在这个示例中，我们将使用二次Lagrange多项式作为基函数，并采用自适应网格细化策略。通过在初始网格上求解局部误差指标，我们可以自动选择需要细化的区域，并将其细化。这样，我们可以在保持高精度的同时有效地减少计算量。
\subsection{代码}
\begin{lstlisting}[language=C++]
// 设置源项函数和边界值函数
class RightHandSide : public Function<dim>
{
public:
  RightHandSide() : Function<dim>() {}

  virtual double value(const Point<dim> &p,
                       const unsigned int component = 0) const override
  {
    return 1.0;
  }
};

class BoundaryValues : public Function<dim>
{
public:
  BoundaryValues() : Function<dim>() {}

  virtual double value(const Point<dim> &p,
                       const unsigned int component = 0) const override
  {
    return p[0]*(1-p[0])*p[1]*(1-p[1]);
  }
};

// 在主函数中进行求解
int main()
{
  //...

  // 定义ProblemDescription对象设置问题参数
  ProblemDescription<dim> problem_description;

  // 设置源项函数和边界值函数
  RightHandSide right_hand_side;
  BoundaryValues boundary_values;
  problem_description.set_rhs_function(right_hand_side);
  problem_description.set_boundary_values(boundary_values);

  // 定义求解器并运行求解
  PoissonSolver<dim> solver(problem_description);
  solver.run();

  //...

  return 0;
}
\end{lstlisting}
\subsection{图片}
在运行step-3程序后，我们可以得到以下计算结果的图片：
\begin{figure}
  \centering
  \includegraphics[scale=0.5]{图片.pdf}
  \caption{PDF的示例图}
  \label{fig:pdf_example}
\end{figure}
以上是对deal.II库示例程序step-3的简要描述。通过离散化Poisson方程并使用有限元方法进行数值求解，我们可以获得边界条件下的数值解。代码示例中的参数设置和实施方法可以用于更一般的问题，并且deal.II库提供了丰富的工具和函数来支持更复杂的数值计算任务。

\end{document}
